Charge fluctuations close to phase separation in the two dimensional t — J model 
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We have studied the t-J model using the Green Function Monte Carlo technique. We have obtained 
accurate energies well converged in the thermodynamic limit, by performing simulations up to 242 
lattice sites. By studying the energy as a function of hole doping we conclude that there is no 
phase separation in the physical region, relevant for HTc superconductors. This finding is further 
supported by the hole-hole correlation function calculation. Remarkably, by approaching the phase 
separation instability, for J c /t ~ 0.5, this function displays enhanced fluctuations at incommensurate 
wavevectors, scaling linearly with the doping, in agreement with experimental findings. 
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One of the most important questions raised in the con- 
text of High-Tc superconductivity is whether the strong 
competition between hole propagation and antiferromag- 
netic (AF) order in the CuOi planes, leads to the segre- 
gation of holes in regions without AF order (phase sep- 
aration) [Q or, analogously, to the formation of one di- 
mensional "stripes" with high hole-density 0. 

The simplest model taking into account spin interac- 
tions and hole kinematics is the t — J model H 

H = J^2(Si ■ Sj - -runs) - t ^ (4c c i° + h.c), (1) 



where cj creates an electron of spin a on the site i, rii 
and S% being the electron number and spin operators re- 
spectively. Double occupations are not allowed and sum- 
mations are extended to nearest neighbors. 

Following [BJ, whenever the ground state energy per 
hole eh (8) as a function of the doping 8 



e h (S) = (e(S)-e )/S 



(2) 



has a minimum for a non zero value of 6 — x c phase 
separation (PS) is energetically stabilized. In fact, by 
using the Maxwell construction, a macroscopic gain in 
energy can be obtained by phase separating the holes 
at all doping 8 < x c in a hole rich phase with 8 = x c 
and a fully AF region without holes with energy eo = 
-1.16944(4) J ffl. 

Many authors |4|j7|-pT|] have tried to clarify the impor- 
tant issue of PS by using numerical techniques. Calcula- 
tions of exact ground state energies on small lattice sizes 
do confirm the existence of PS at physical J/t and 
doping , but the sizes considered are far to be represen- 
tative of the thermodynamic limit, where the question 
is meaningful. By contrast, using the high temperature 
expansion (HTE), PS was found only for large J/t 11 1. 

Much more information, not only on PS, can be ob- 
tained by studying the behavior of the hole - hole corre- 
lation function: 



N(q) 



where L is the number of lattice sites. The onset of PS 
is characterized by the divergence of N(q) for small mo- 
menta q. According to the HTE (l2), N(q) displays some 
interesting feature at twice the spinless fermion Fermi 
vector but no PS in the physical region. These results 
are not affected by small size effects but are obviously 
limited by the difficulty of extrapolating an high temper- 
ature series to zero temperature. 

In this letter we tackle the PS problem by the lattice 
Green Function Monte Carlo (GFMC) technique. This 
method allows to filter out from a given trial wavefunc- 
tion \i(jt) the ground state \if>o) of H by statistical appli- 
cation of the power method: 



^o)oc(A-iJ)"|^ 



(4) 
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where A is a suitably large constant required to allow con- 
vergence to the ground state for large n. In this scheme 
a Markov process is defined, acting on electron config- 
urations {x} with definite positions and spins, which 
are changed according to the Green Function G x ' iX = 
A£ x ',a; — H x > lX [^3). The algorithm is efficient even for 
large system size, provided all the G x ', x are non negative, 
otherwise one is facing the well known " sign problem" . 
For the diagonal elements the above condition can be 
easily fulfilled by increasing A. However an exceedingly 
large value of A -which is often the case especially for 
fermions- determines a slowing down of the algorithm, 
since there is a very small probability ~ 1/ A to accept a 
new configuration x' and the algorithm remains almost 
always stacked in the old one x. Thus one needs much 
more power iterations (|j) to generate statistically inde- 
pendent configurations. In order to overcome this dif- 
ficulty, following JhJ, it is better to determine a priori 
the number of diagonal moves before an off diagonal is 
accepted. Thus one can generate each time a new con- 
figuration without caring of a very large value of A. A 
further improvement is to let A — > 00 with n — Ar, so 
that Eq. (^J) becomes a more conventional imaginary time 
evolution by the operator e~ Hr . Further details of this 
technique will be given elsewhere. 

The negative off-diagonal elements of the Green Func- 
tion can be handled approximately by a recent achieve- 
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ment in GFMC: the fixed node technique (FN) for lat- 
tice Hamiltonians [ p"5| |. In order to control and improve 
the accuracy of the FN approximation in a systematic 
way, we have used a recent development, the GFMC with 
stochastic reconfiguration (GFMCSR) @. The GFM- 
CSR results presented here have been obtained by keep- 
ing unchanged before and after each stochastic recon- 
figuration the energy, the kinetic energy and the nearest 
neighbor hole-hole correlation function. GFMCSR has to 
be implemented with a very large number M of walkers 
in order to stabilize the simulation (M ~ 2000 -r 4000), 
whereas in the FN calculations M — 800 is by far suffi- 
cient to control the bias due to the finite walker popula- 
tion with the technique described in Ref. ||. 
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FIG. 1. eh{5) computed with VMC (long dashed lines), 
FN (short dashed lines), GFMCSR (dotted lines) compared 
with the exact results |20| (full lines) for a 26 lattice size. 
Error bars are much smaller than the size of the symbols. 
Lines (guides to the eye) connect the two and four hole results. 

A fundamental ingredient in GFMC is the choice of 
the guiding function in order to perform importance sam- 
pling, as described in ||. At finite doping we have used 
a pure d-wave BCS guiding wave function pTq— [lq] plus 
a long range density-density Jastrow factor p9| . At half 
filling instead we use the guiding function described in || 
that allows to obtain the exact answer for eo, as there is 
no sign problem at zero doping for this particular guiding 
function. In the following we use the exact eo at the given 
size L for all variational (VMC), FN and GFMCSR cal- 
culations of eu (S). Further we employ periodic boundary 
conditions tilted by 45 (~ 11) degrees on square lattices 
with L = 2l 2 and odd integer I > 5 (L = 26). 

We have compared our Monte Carlo calculations with 
the exact Lanczos results for the largest size L = 26 
available in the literature |^(J. We have found that the 
FN approximation improves the ground state energy of 
the best starting variational (and guiding) wavefunction 
by a factor of three (Fig. g) and the GFMCSR by an- 
other similar factor, yielding finally an accuracy of less 
than 100K on the energy per hole, which is physically 
acceptable if compared with the low energy coupling of 
the model J ~ 1500iV 
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FIG. 2. eh(S) computed with FN (dashed lines) and 
with GFMCSR (continuous lines) for L = 50 (triangles) and 
L = 98 (circles) at J = 0.4t. Lines are guides to the eye. 

As shown in Fig. [l] this kind of accuracy depends 
weekly on the number of holes. However, for small lat- 
tices, the main difficulty to detect PS is the resolution 
in doping. By increasing the system size (see Fig. ||), 
the difference between FN calculation and the GFMCSR 
one remains of the same order, and much below the VMC 
energies. Thus we expect that the accuracy of the cal- 
culation is not very much size dependent, even for large 
systems where no exact solution is available. We remark 
that, as the accuracy of the calculation is improved from 
FN to GFMCSR, the minimum in the hole energy dis- 
appears for the largest size in Fig. ||. This suggests that 
the occurrence of PS at J/t = 0.4 and at this system size 
is an artifact of the FN. This approximation naturally 
enhances the tendency of PS because the FN acts only 
on the kinetic part of the Hamiltonian, thus implying a 
tendency to localize the holes. Nonetheless, as shown in 
Fig. |^, we have obtained, even within the FN, the stabil- 
ity of the uniform phase in the thermodynamic limit at 
J/t = 0.4. Thus we conclude that PS at this J/t value is 
only a finite size effect. 
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FIG. 3. eh(S) computed with FN for several lattice sizes. 

Several authors |7]|| have used the infinite L limit for 
eo even for the finite L evaluation of e/j (S) . We have in- 
stead used the exact eo for each lattice size, and checked 
that for the largest size calculations both choices of eo 
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lead to the same e/ 1 (<5). Indeed the small lattice size re- 
sults are very sensitive at low doping to the particular 
choice of eo, and this may explain the contradictory re- 
sults presented in the literature so far. Moreover in Ref. 
H the most important and delicate low doping region is 
studied only with fairly small lattice sizes. 

The absence of PS is also confirmed by the low mo- 
menta behavior of N(q), which will be discussed in the 
remaining part of this letter. In order to compute the 
ground state correlation functions N(q) we have used 
two different methods: the "forward walking technique" 
(FW) which allows the direct evaluation of the ground 
state expectation value, at the expense of very large er- 
ror bars when the convergence to the ground state is 
particularly slow, as shown in (H); the second technique 
is based on Hellmann-Fcynman theorem, and amounts to 
compute the ground state energy E(A) of the Hamilto- 
nian in presence of a small perturbation XN(q), N(q) = 
-jjE(A) |a=0i the derivative being estimated numerically 
by a few runs for different A's. The latter technique is 
much more stable, especially for large size, but each q 
value requires several simulations, whereas a single one 
is sufficient for the FW technique for all g's. Thus we 
have used the more expensive method for the small q val- 
ues where the FW convergence is more difficult, and we 
have checked the consistency of both methods in the re- 
maining momentum region. Moreover the FW technique 
is limited to the FN approximation and has not been 
extended within the GFMCSR scheme yet. For the 26 
site cluster the FN results for N(q) are accurate within 
3%, as compared with the exact diagonalization data. 
Since this accuracy is already very good for determin- 
ing the qualitative features of N(q), we have chosen to 
work within FN, avoiding to implement the much more 
expensive GFMCSR on large clusters. 

In Fig. f|the N(q) is plotted for several dopings and lat- 
tice sizes. Inside the PS region, at J = t (Fig. [| b), N(q) 
shows a divergent peak at small q, as can be expected 
from general arguments. Instead for J = OAt (Fig. [| a) 
the charge correlations approach zero as q — > 0, confirm- 
ing the absence of PS even at the lowest doping. Nev- 
ertheless enhanced fluctuations are clearly evident along 
the (1, 1) and (1,0) directions. We believe that these in- 
commensurate peaks are a genuine feature of the ground 
state of the model, because they do not appear for in- 
stance at the VMC level, and it is extremely important 
to use many power iterations (^) to eliminate the bias due 
to the VMC guiding wavefunction. The N(q) at this J/t 
value is very weakly size dependent, so that the overall 
shape of this function in the thermodynamic limit should 
not differ too much from the one shown in Fig. || (a). 
Thus N(q) should be always finite even for small q, ruling 
out PS and charge density wave instability. Even though 
we found some peaks at incommensurate wavevectors, 
that maybe reminiscent of some dynamical stripe order, 
they are not consistent with a static stripe structure (the 



peaks should diverge). 
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FIG. 4. (a) N(q) for several doping 8 — 0.06 (triangles), 
8 — 0.10 (squares), 8 = 0.17 (circles) for different sizes, L — 98 
(full symbols) and L — 162 (empty symbols). Incommensu- 
rate peaks are shown by the arrows. V — (0,0), X — 
M = (-7T, 0). Inset: empty dots (full dots) represent peak posi- 
tions along the (1,0) direction in q space (diagonal direction). 
The dashed line displays the 4tt slope, (b) N(q) for 8 — 0.06 
and L = 162 (empty dots), L = 98 (full dots) 

In fact we tried to stabilize some static stripe order 
similarly to what was done within Density Matrix Rcnor- 
malization Group (DMRG) |^l| , but in our approach the 
uniform phase remains always the most stable one at 
J/t = 0.4. The stabilization of the stripe phase with 
DMRG is probably due to the use of open boundary 
condition in one direction. We have not attempted to 
use open boundary condition as, i) the momentum is not 
defined on finite lattices, ii) it is more difficult to reach 
the thermodynamic limit, iii) the gap to the first excited 
state is much reduced, especially close to a PS instability, 
yielding a much slower convergence of the power method 

Remarkably, as S is increased the peak at finite mo- 
mentum moves far from the T point at a distance that 
scales linearly with the doping with a coefficient which is 
surprisingly close to 4ir. This is exactly the coefficient ob- 
tained experimentally in La2- x -sNd x SrsCuOi |22|. We 



have also found that this peak position does not depend 
on J/t, implying that the 4ir5 slope could be a general 
feature of the t — J model in the region without PS. It 
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is reasonable to expect that the interaction of electrons 
with the ions could further enhance the intensity of this 
incommensurate peaks, leading not only to a qualitative 
but also to a quantitative agreement with experiments 

ill- 

In conclusion we have performed accurate calculations 
of the ground state energy of the 2D t — J model with 
three methods: VMC, FN and GFMCSR. 



and ii) long range forces to push the q — PS instability 
to finite q value Q . 
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FIG. 5. The low doping PS phase diagram of the t — J 
model. Full dots (connected by continuous line) are computed 
with FN, the full triangle with GFMCSR, empty dots are from 
ji| and the dashed line is from 

As summarized in the phase diagram picture (Fig. |^) 
we find no evidence of PS for J < OAt, at least for 
5 > 2.5%, which represents the smallest doping consid- 
ered here (6 holes in a L = 242 sites), and a transition to 
the phase separated regime at a critical J c ~ 0.5i. The 
small doping difference between FN and GFMCSR for 
the occurrence of PS at J = t strengthens the validity 
of the FN phase separation diagram. These results are 
in acceptable quantitative agreement with Ref. RPlJlOl 
but we believe that our calculation represents a much 
better attempt to control the finite size effects, which 
are very important, especially at small doping. Instead, 
in the large doping region, the best agreement is found 
with the exact diagonalization data Q] on small clusters. 
Close to PS we have computed the equal time charge cor- 
relations and enhanced fluctuations at incommensurate 
momenta have been found. The fact that the position of 
these peaks approach the T point linearly with the dop- 
ing cannot be explained within a simple model for the 
holes, like e.g. the spinless fermion one, as was proposed 
by |Q to explain the shape of N(q) using HTE. In this 
case in fact the characteristic incommensurate wavevec- 
tor 2kp approaches the T point in a much more singular 
fashion 2kp ~ \f8. On the other hand the hard-core bo- 
son model is unable to produce any feature at momenta 
different from T. Further study is probably necessary 
to clarify this point. We emphasize however that, within 
the t — J model, it is possible to reproduce the qualitative 
features of the experimental findings of incommensurate 
charge fluctuations, without i) electron-phonon coupling 
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